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5 Serial No. 60/1 5 1 ,61 8, filed 08/3 1/1 999, the entire text of which is hereby incorporated by 
reference. 
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Nos. ACI-9624957, ACI-972 1 349, ASC-89202 1 9 awarded by the National Science 
Q 10 Foundation; Grant No. DMS-9875042 awarded to Caltech's OPAAL Project by DARPA and 
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H TECHNICAL FIELD 

IJI This invention relates to computer techniques for performing finite element analysis, 

U 

and more particularly to novel computer techniques for performing finite element analysis of 
15 thin-shell models using subdivision surfaces. 



m 



BACKGROUND 

The Kirchhoff theory of thin plates and the Kirchhoff-Love theory of thin shells are 
characterized by energy fimctionals w^hich depend on curvature; consequently they contain 
second-order derivatives of displacement. The resulting Euler-Lagrange - or equilibrium - 

20 equations in tum take the form of fourth-order partial differential equations. It is well-known 
from approximation theory that in this context, the convergence of finite-element solutions 
requires so-called interpolation. More precisely, in order to ensure that the bending energy 
is finite, the test functions have to be if, or square-integrable functions whose first- and 
second-order derivatives are themselves square-integrable. Unfortunately, for general 

25 unstructured meshes it is not possible to ensure continuity in the conventional sense of 
strict slope continuity across finite elements when the elements are endowed with purely 
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local polynomial shape functions and the nodal degrees of freedom consist of displacements 
and slopes only. Inclusion of higher-order derivatives among the nodal variables leads to 
well-known difficulties, e.g., the inability to account for stress and strain discontinuities in 
shells whose properties vary discontinuously across element boundaries, and, owing to the 
5 high order of the polynomial interpolation required, the presence of spurious oscillations in 
the solution. 

The difficulties inherent in interpolation have motivated a number of altemative 
approaches, all of which endeavor to "beat" the continuity requirement. Examples are: 
quasi-conforming elements obtained by relaxing the strict Kirchhoff constraint; the use of 
10 Reissner-Mindlin theories for thick plates and shells (which requires conventional 

interpolation only); reduced-integration penalty methods; mixed formulations; degenerate 
solid elements; and many others known from the literature. C° elements often exhibit poor 
performance in the thin-shell limit - especially in the presence of severe element distortion. 
^ Such distortion may be due to a variety of pathologies such as shear and membrane locking. , 

|=i 15 The proliferation of approaches and the rapid growth of the specialized literature attest to the 



m 



inherent, perhaps insurmountable, difficulties in vanquishing the C continuity requirement 



p SUMMARY 

^= The invention includes a method, system, and computer program for performing finite 

element analysis on a shell based on the use of subdivision surfaces for: (1) describing the 

Q 

Q 20 geometry of a shell in its undeformed configuration, and (2) generating smooth interpolated 
displacement fields possessing bounded energy within the strict framework of the Kirchhoflf- 
Love theory of thin shells. The preferred subdivision strategy is Loop's scheme, with 
extensions to account for creases and displacement boundary conditions. The displacement 
fields obtained by subdivision are If and, consequently, have a finite Kirchhoff-Love energy. 
25 The resulting finite elements contain three nodes and element integrals are computed by a 
one-point quadrature. The displacement field of the shell is interpolated from nodal 
displacements only. In particular, no nodal rotations are used in the interpolation. The 
interpolation scheme induced by subdivision is nonlocal, i.e., the displacement field over one 
element depends on the nodal displacements of the element nodes and all nodes of 
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immediately neighboring elements. However, the use of subdivision surfaces ensures that all 
the local displacement fields thus constructed combine conformingly to define one single 
limit surface. Numerical tests, including known obstacle courses of benchmark problems, 
demonstrate the high accuracy and optimal convergence of the method. 

More particular, in one aspect, the invention includes a method, system, and computer 
program for performing finite element analysis on a shell by the following steps: 

• Modeling the geometry of the shell using subdivision surfaces. 

• Characterizing an environment for the shell, including environmental factors affecting 
the mechanical behavior of the modeled shell (such environment factors includes 
loading conditions, material properties, and boundary conditions for the modeled 
shell; the loading conditions may include an indication of applied forces and of 
thermal loading). 

• Computing the mechanical response of the modeled shell, taking into account the 
characterized environment, using a finite element analysis; the finite element analysis 
preferably uses subdivision basis functions as shape functions, but may in general use 
any suitably smooth shape function. 

• Outputting a description of the geometry of the modeled shell as determined from the 
computed mechanical response (optionally, the output can include outputting 
indications of the characterized environment). 

A specific embodiment of the inventive method for performing finite element analysis 
using subdivision basis functions can be summarized as follows: 

• input a mesh comprising a set of data points each having connectivity to neighboring 
data points, the mesh defining physical parameters; 

• specify an initial state for the mesh; 

• define a set of linear differential equations comprising a stiffness matrix and an 
external forcing vector, at least one such equation having a fourth order differential 
operator; 

• solve the set of linear equations as applied to the mesh; 
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• output the solution to the set of linear equations as defining a modification of the 
initial state of the mesh based on the stifftiess matrix and in response to the external 
forcing vector. 

This procedure can be implemented by programming the relevant equations set forth 
below in light of known techniques for finite element analysis. 

The details of one or more embodiments of the invention are set forth in the accompa- 
nying drawings and the description below. Other features, objects, and advantages of the 
invention will be apparent fi^om the description and drawings, and from the claims. 

DESCRIPTION OF DRAWINGS 

FIG. 1 is a diagram showing a shell geometry in a reference configuration and a 
deformed configuration. 

FIG. 2 is a diagram showing graphically depicting refinement of a coarse triangular 
mesh to a finer mesh by quadrisection. 

FIG. 3a is a diagram showing a simple one-dimensional example of an approximating 
subdivision scheme. FIG. 3b is a diagram showing an interpolating subdivision scheme. 

FIG. 4a is a diagram showing a subdivision mask for computing x\ . FIG. 4b is a 
diagram showing a subdivision mask for computing Xj , ... . 

FIG. 5 is a diagram showing a graphical example of a boundary approximation using 
temporary vertices. 

FIG. 6 is a diagram showing an example of a geometrically complex shape (a 
distributor cap), demonstrating the ability of subdivision surfaces to model intricate shapes in 
a straightforward fashion. 

FIG. 7a is a diagram showing a limit position mask. FIGS. 7b and 7c are diagrams 
showing two tangent masks corresponding to the limit position mask shown in FIG. 7a. 

FIG. 8 is a graphical example of a regular box-spline patch v^th 12 control points. 

FIG. 9 is a diagram showing the basic idea for fianction evaluation at arbitrary 
parameter values. 
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FIG. 10 is a diagram showing the relation between the coordinates (9*, 6") of original 
triangles and the coordinates (9\ 6^) of refined triangles. 

FIGS. 11a, lib, and 1 Ic are diagrams showing various boundary conditions for the 
cases of displacements and rotations fixed (FIG. 1 la), displacement and rotations fi-ee (FIG. 
1 lb), and displacements fixed and rotations fi-ee (FIG. 11c). 

FIG. 12 is a flowchart showing a functional method for performing finite element 
analysis using subdivision basis functions. 

FIG. 13 is a flowchart showing one analytical solution process for performing finite 
element analysis using subdivision basis fiinctions. 

FIG. 14 is a diagram shov^dng an example of a rectangular plate. 

FIGS. 15a and 15b are diagrams respectively showing the computed limit surfaces of 
the simply-supported plate of FIG. 14 and the clamped plate of FIG. 14 following 
deformation. 

FIGS. 16a and 16b are graphs respectively showing the convergence of the computed 
maximum-displacement and energy-norm errors as a fimction of the number of degrees of 
freedom for the plate of FIG. 15b. 

FIG. 17 is a diagram showing an example of a Scordelis-Lo Roof cylindrical shell. 

FIGS. 18a and 18b are diagrams respectively showing the undeformed control mesh 
and the deformed limit surface for the shell of FIG. 17. 

FIG. 19 shows a convergence plot for the maximum displacement of the deformed 
limit surface of FIG. 18b. 

FIG. 20a is a diagram showing an example of a cylindrical shell. FIG. 20b is a 
diagram showing a typical control mesh for the cylinder of FIG. 20a. FIG. 20c is a diagram 
showing the computed deformed limit surface for the cylinder of FIG. 20a, with two point 
loads . 

FIGS. 21a and 21b are graphs respecfively showing displacement-convergence results 
for the rigid diaphragm case and the fi-ee-end case for the cylinder of FIG. 20a, FIG. 20b, and 
FIG. 20c. 
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FIG. 22a is a diagram showing an example of a hemispherical shell. FIGS. 22b and 
22c are diagrams respectively showing a typical control mesh and the corresponding 
deformed limit surface for the hemispherical shell of FIG. 22a. 

FIG. 23 shows the convergence of the radial normalized displacement under the 
5 applied loads for the hemispherical shell of FIG. 22a, FIG. 22b, and FIG. 22c. 

Like reference numbers and designations in the various drawings indicate like 
elements. 



DETAILED DESCRIPTION 

10 Introduction 

Simultaneously with the development of plate and shell elements, and for the most 
part unbeknownst to mechanicians, the field of computer aided geometric design has taken 
considerable strides towards the efficient generation and representation of smooth surfaces. 
In particular, the use of subdivision surfaces provides a powerful tool for generating smooth 

15 surfaces which either interpolate or approximate an arbitrary collection of points or "nodes". 
Here, smoothness is understood in the sense of }f surfaces, /.e., surfaces whose curvature 
tensor is Zj, or square summable. Subdivision surfaces follow as the limit of a recursive 
refinement of the triangulation of the nodal point set, e.g., by recourse to the classical Loop 
scheme. Within this framework, the treatment of complex geometries with intersections or 

20 curved boundaries is straightforward. Subdivision surfaces obtained by the Loop scheme are 
guaranteed to be if, i.e., to have finite bending energy, and are therefore ideally suited as test 
functions for plate and shell analyses. The smoothness of the limit surface may also be 
suitably relaxed in the presence of thickness or material discontinuities. The method of 
subdivision surfaces thus effectively solves the elusive and long-standing -interpolation 

25 problem which has traditionally plagued plate and shell finite-element analyses. The ready 
availability of smooth approximating surfaces for arbitrary geometries and triangulations 
enables a return to the "basic" finite-element method, e.g., the Rayleigh-Ritz method, with 
the attendant guarantee of optimal convergence. 
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We teach the subdivision-surface concept as a new paradigm for plate and shell C* 
finite-element analyses. We rely on subdivision to generate smooth deformed surfaces from a 
triangulation of an arbitrary nodal point set. This nodal point set is displaced from that which 
defines the reference configuration of the shell according to an array of nodal displacements. 

5 The nodal displacements are determined as follows. The energy of the deformed plate or shell 
is given by a direct evaluation of the Kirchhoff or Kirchhoff-Love energy fimctional. The 
requisite boundedness of the bending energy is ensured by the H~ property of the test 
deformed geometries. The equilibriimi displacements then follow simply by recourse to 
energy minimization. Within the framework of linear theories, this process of energy 

10 minimization leads to a synmietric and banded system of linear equations for the nodal 
displacements. 

The triangles in the triangulation of the nodal point set may be regarded as three-node 
finite elements. In particular, the total energy of the shell is the sum of the local energies of 
the elements. These local energies in tum follow by integration over the domain of the 

15 element. However, the interpolation scheme to which the subdivision paradigm leads differs 
from conventional finite-element interpolation in a crucial respect: the displacement field 
within an element depends not only on the displacements of the nodes attached to the element 
but also on the displacements of all the immediately adjacent nodes in the triangulation. 
Thus, the displacement field wdthin an element is determined by the nodal displacements of a 

20 "link" or "patch" of adjacent elements. (However, special rules are required for elements 

which abut on an edge of the shell.) Our approach shares some aspects in common with the 
finite-volume approach recently proposed by J. Rojek, E. Onate, and E. Postek, "Application 
of explicit fe codes to simulation of sheet and bulk metal forming processes," J. Mater. 
Process. Techn., 80:620-627, 1998. For instance, the patches corresponding to neighboring 

25 elements may overlap. However, when the displacement field is constructed by subdivision, 
as in our approach, the displacement representations within the intersection of two patches 
coincide exactly. Thus, subdivision leads to a unique and well-defined surface over the 
complete finite-element triangulation, as opposed to a collection of non-conforming local 
interpolations. 
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An additional advantage afforded by the present approach is that the geometrical 
modeling and the finite-element analysis are based on an identical representational paradigm, 
that is, both the undeformed and the deformed geometries of the plates and shells are 
described by recourse to subdivision. The prior art has stressed the necessity of a unique 
5 framework for both geometric design and mechanical analysis. In contrast, the unification of 
the geometrical and finite-element representations in accordance with the present invention 
offers a robust environment which effectively sidesteps many of the difficulties inherent in 
the currently available software tools, which suffer fi-om heterogeneous and, therefore, error- 
prone interfaces. 

10 Concisely stated, the invention includes a method, system, and computer program for 

performing finite element analysis on a shell by the following steps: 
(a) Modeling the geometry of the shell using subdivision surfaces. 
^ (b) Characterizing an environment for the shell, including environmental factors 

=P affecting the mechanical behavior of the modeled shell (such environment factors 

15 includes loading conditions, material properties, and boundary conditions for the 

H . 5 

modeled shell; the loading conditions may include an indication of applied forces as 
well as thermal loading). 

(c) Computing the mechanical response of the modeled shell, taking into account the 
characterized environment, using a finite element analysis; the finite element 

20 analysis preferably uses subdivision basis functions as shape functions, but may in 

general use any suitably smooth shape function. 

(d) Outputting a description of the geometry of the modeled shell as determined from 
the computed mechanical response (optionally, the output can include outputting 
indications of the characterized environment). 

25 The Thin-Shell Boundary Value Problem 

In this section we summarize the field equations for the classical stress-resultant shell 
model. Here we follow the elegant formulation by Simo et al (J.C. Simo and D.D. Fox, "On 
a stress resultant geometrically exact shell model, part i: Formulation and optimal 
parameterization," CompuL Methods Appi Meek Engrg,, 72:267-304, 1989; J.C! Simo, D.D. 
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Fox, and M.S. Rifai, "On a stress resultant geometrically exact shell model, part ii: The linear 
theory; computational aspects," Comput, Methods Appl Mech, Engrg., 73:53-92, 1989) of 
the Reissner-Mindlin theory, which we specialize to Kirchhoff-Love theory by explicitly 
constraining the shell director to remain normal to the deformed middle surface of the shell. 
Throughout the present work we confine our attention to the linear theory of shells under 
static loading. In Section we briefly summarize the relevant aspects of the standard finite- 
element discretization of the Kirchhoff-Love thin-shell theory. 

Kinematics of Deformation 

FIG. 1 shows a shell geometry in a reference configuration 12, a deformed 
configuration 14, and the underlying parameter space 10. The undeformed geometry 12 is 
characterized by a middle surface of domain □ and boundary ID = 50. For simplicity we 
assume throughout that the thickness h of the shell is uniform. The undeformed shell 12 
deforms under the action of applied loads (which may include thermal loading) and adopts a 
deformed configuration 14, characterized by a middle surface of domain Q and boundary Y = 

da. 

The position vectors r and r of a material point in the reference 12 and deformed 
configurations 14 of the shell may be parameterized in terms of a system {Q\ 9^,0^} of 
curvilinear coordinates as: 

The functions jc (9\9^) and jc(9',9^) furnish a parametric representation of the middle surface 
in the reference 12 configuration and deformed configuration 14. The corresponding surface 
basis vectors are: 
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where here and henceforth Greek indices take the values 1 and 2, and a comma is used to 
denote partial differentiation. The covariant components of the surface metric tensors in turn 
follow as: 

whereas the covariant components of curvature tensors are: 

For later reference we also introduce the contravariant components of the undeformed 

and deformed surface metric tensors, a"^ and a^^ , respectively. The defining property of 
these components is: 

We also note that the element of area over □ follows as 
where 

Va = \ai X 

is the jacobian of the surface coordinates {9*, 6^}. 

The shell director as in the reference configuration 12 coincides with the normal to 
the undeformed middle surface of the shell and hence has the properties: 

• ^3 = 0 J \as\ = 1 
which give as explicitly in the form: 
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as = Tzz (10) 

\ai X a^l 

For the moment we allow the director on the deformed configuration 14 of the shell to be 
an arbitrary vector field. 

The covariant base vectors in the reference and the current configurations 12, 14 
follow simply as: 

^a— ^ = + ^^«3,a , ^3 = ^ = ^3 (H) 

fl^a = ^ = «a + ^^«3,a , = = (12) 

The corresponding covariant components of the metric tensors in both configurations are: 

9ij = 9i ' 9j , 9ij = 9i • 9j (13) 

where here and henceforth lowercase Latin indices take the values 1, 2, and 3. 

The Green-Lagrange strain tensor is defined as the difference between the metric 
tensors on the deformed and undeformed configurations of the shell, i.e., 

Eij^^{9ij~gij) (14) 
Using Eqs. 11, 12, and 13, the Green-Lagrange strain of the shell is found to be of the 

form: 

Eij = aij+9^d.,, (15) 

to first order in the shell thickness h. The nonzero components of the tensors a,-, and P,-^ are in 
turn related to the deformation of the shell as follows: 

^ij = \i^i' ~ * ^j) (16) 
PaP = CLa • a3./3 " * a3,/3 ' (17) 

In particular, the in-plane components a^p, or membrane strains, measure the straining of the 
surface; the components a^^ measure the shearing of the director as ; the component 
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measures the stretching of the director; and the components P^.^, or bending strains, measure 
the bending or change in curvature of the shell, respectively. 

The above kinematic relations allow for finite deformations as well as for shearing 
and stretching of the shell director. In the remainder of this description, we restrict our 
attention to the Kirchhofif-Love theory of thin shells and, accordingly, we constrain the 
deformed director to coincide with the unit normal to the deformed middle surface of the 
shell, /.e., 

• as ^ Oj \az\ = 1 
which yields: 

ai X a2 



as = 



\ai X a2| 



Owing to these constraints, the shear strains a^,3 vanish identically and the bending 
strains simplify to: 

It follows from these relations that, by virtue of the assumed Kirchhoff-Love kinematics, all 
the strain measures of interest may be deduced from the deformation of the middle surface of 
the shell. 

For simplicity, we restrict the scope of subsequent discussions to linearized 
kinematics. To this end, we begin by v^iting 



where u(Q\ 0^) is the displacement field of the middle surface of the shell. To first order in u 
the membrane and bending strains then follow as: 
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1/- - X 



+ 1^,1 -(02 X as) + n,2 -(aa x ai)J 



It is clear from these expressions that the displacement field u of the middle surface furnishes 
a complete description of the deformation of the shell and may therefore be regarded as the 
primary unknown of the analysis. It also follows that the deformed and undeformed domains 
Q and O are indistinguishable to within the order of approximation of the linearized theory 
and, in consequence, we drop the distinction between the two domains throughout the 
remainder of this description. 

Equilibrium Deformations of Elastic Shells 

Next, we seek to characterize the equilibrium configurations of the shell by recourse 
to energy principles. For simplicity, we shall assume throughout that the shell is linear elastic 
with a strain energy density per unit area of the form: 

^(«' /3) = I i^^'^^^^^.^a,, + l^^0-^^H'^^^0^,P,s (24) 

where E is Young's modulus, v is Poisson's ratio, and 

= i/a^V* + i(l - u) {aTaP^ + al'^ot') (25) 

In Eq. 24, the first term is the membrane strain energy density and the second term is the 
bending strain energy density. The membrane and bending stresses follow from Eq. 24 by 
work conjugacy, with the result: 

= ^ = T^^H'^'^'cc^s (26) 

The membrane and bending stress tensors rf-^ and may be given a direct mechanistic 
interpretation as force and moment resultants. 
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The shell is subject to a system of external dead loads consisting of distributed loads q 
per unit area of Q and axial forces A'^ per unit length of P. Under these conditions the 
potential energy of the shell takes the form: 

$[u]-$^"'[n] + $*^[ti] (28) 

where 

^^^[u] = [ W{ct.p)dn (29) 

is the elastic potential energy and 

^ - f q-udQ- f N ^uds (30) 

is the potential energy of the applied loads. 

The stable equilibrium configurations of the shell now follow from the principle of 
minimum potential energy: 



m $M=mf^$M (31) 

j; — z 

where Fis the space of solutions consisting of all trial displacement fields v with finite 
ry 10 energy 0[v]. It is clear from the form of the elastic energy of the shell that such trial 
e displacement fields must necessarily have square integrable first and second derivatives. 

''^ Within the context of the linear theory, and under suitable technical restrictions on the 

domain Q and the applied loads, it therefore follows that Fmay be identified with the 
Sobolev space of functions H\Q, R^). In particular, an acceptable finite-element interpolation 
15 method must guarantee that all trial finite element interpolants belong to this space. 

The Euler-Lagrange equations corresponding to the minimum principle (Eq. 31) may 
be expressed in weak form as:. 



-14- 



Docket No.: 0661 8/5QiMl/CIT3061 



{D^u],Su} = {D^'"'[u],du) + {D^'''^[u]Ju) = 0 

which is a statement of the principle of virtual work. Here {D 0[u],du > denotes the first 
variation of O at a in the direction of the virtual displacements du. 



is the internal virtual work and 



{D^^^[ul6u) = - j^q ' 5udQ - j^N ' duds 



is the external virtual v^ork. The minimum principle (Eq. 3 1) or, equivalently, the virtual 
work principle (Eq. 32) are subsequently taken as a basis for formulating finite-element 
approximations to the equilibrium configuration of the shell. 

Finite-Element Discretization 

We now turn to the finite-element discretization of the potential energy of the shell 
(Eq. 28) or, equivalently, its first variation (Eq. 32). To this end, it proves convenient to 
adopt Voigt's notation and map symmetric second-order tensors into arrays by recourse to 
the conventions: 

^11 ^ ^11 >^ / au \ ( 0n 




(322 



The constitutive relations, Eqs. 26 and 27, may similarly be written in the form: 
Eh 

n = - -Hot 

1 -- 



where we write: 
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if = 



/ (a")2 2/aiia22 + (1 - z/)(ai2)2 
\ sym. 



which replaces Eq. 25 within the Voigt formalism. Using these conventions, the internal 
virtual work (Eq. 33) may be recast in the convenient form 



(^■>],(Jn)= / 

Jf2 



Eh 



Eh!" 



12(1 - i/2) 



Next we proceed to partition the domain Q of the shell into a finite element mesh, the 
precise nature of which will remain unspecified for now. The collection of element domains 
in the mesh is {Q^^, K = 1 ,...J^EL}, where denotes the domain of element K and NEL is 
the total number of elements in the mesh. The finite-element mesh may be taken as a basis 
for introducing a displacement interpolation of the general form: 

1=1 

where {V, /= 1,...,A^ } are the shape functions, {«/,/= 1,...,A/^ } are the corresponding nodal 
displacements, and A'^ is the number of nodes in the mesh. Owing to the linearity of the 
dependence of the finite-element interpolant m^, on the nodal displacements, an application of 
Eq. 22 and Eqs. 23 to 40 gives the finite-element membrane and bending strains in the form: 

1=1 

= f:B'{9\e')uj 

1=1 

for some matrices and B^. The precise form of these matrices is given in Appendix A.0.3 
below. Finally, the introduction of Eqs. 40, 41, and 42 into the principle of virtual work (Eq. 
32) yields the equations of equilibrium for the nodal displacements: 
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KhUn = /fc (43) 
where, by a slight abuse of notation, we take to signify the array of nodal displacements. 



NEL . 
K=l 



NEL 

1-1/2' ^ 12(1-1/2) 



K=l 



is the stiffiiess matrix, and 

NEL , . . iVEL 

= E L / nnUs} ^ y: fK (45) 

is the nodal force array. It should be carefully noted that, as expected, the global stifftiess and 
force arrays just defined follow by the assembly of low-dimensionality element stiffness and 

force arrays and // , respectively, as is standard in the finite-element method. 

The computation of element arrays requires the evaluation of integrals extended to the 
domain of each element, cf Eqs. 44, 45. These integrals may efficiently be evaluated by 
recourse to numerical quadrature without compromising the order of convergence. For 
instance, the application of a quadrature rule to the calculation of the element stiffness 
matrices leads to. the expression: 



G=l L 



r Eh ; Efl' 



where (G'^, 0^g) the quadrature points, Wq are the corresponding quadrature weights, and 
NQ is the number of quadrature points in the rule. The element force arrays // may be 
computed likewise. 

Sufficient conditions for the quadrature rule to preserve the order of convergence of 
the finite-element method may be found in G. Strang and G. J. Fix, An Analysis of the Finite 
Element Method, Prentice-Hall, Englewood Cliffs, N.J., 1973. In general, this conditions 
demand that certain shape function derivatives be computed exactly and, consequently, place 
a lower boimd on the order NQ of the quadrature rule. These theoretical considerations and 
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our numerical tests show that a one-point quadrature rule is sufficient to compute all element 
arrays of interest. 

Subdivision Surfaces and Finite-Element Interpolation 

It is clear from the preceding developments that the central problem in thin-shell 
finite-element analysis is the formulation of shape functions which are (or "C*" in the 
usual finite-element terminology), as this property ensures the finiteness of the energy of the 
trial displacement fields. In this section we develop one such interpolation scheme based on 
the notion of subdivision surface. 

Subdivision Surfaces 

We begin by reviewing the essential ideas behind subdivision surfaces using one- 
dimensional examples first, then moving to the two-dimensional manifold (with boundary) 
setting relevant to shells. Here we limit ourselves to reviewing various elementary properties 
of subdivision. 

Throughout this discussion the term vertex will be used to refer to nodes in the mesh. 
For example, a vertex has associated with it a particular nodal position. A vertex has a 
topological neighborhood defined by the structure of the mesh. For example, its l-ring 
neighbors are all those vertices which share an edge with it (and recursively for its A:-ring). 
This distinction between vertices and nodal positions is typically not needed when dealing 
with finite elements, but it is important to keep these distinctions in mind when dealing with 
smooth subdivision. 

At the highest level of description we may say that subdivision schemes construct 
smooth surfaces through a limiting procedure of repeated refinement starting from an initial 
mesh. This initial mesh will also be referred to as the control mesh of the surface. Generally, 
subdivision schemes consist of two steps. First the mesh is refined, e.g., all faces are 
quadrisected, followed by the computation of new nodal positions. These positions are 
simple, linear fiinctions of the nodal positions of the coarser mesh. See FIG. 2, which 
graphically depicts refinement of a coarse triangular mesh 20 to a finer mesh 22 by 
quadrisection. For the schemes of interest these computations are local, i.e., they involve only 
nodal positions of the coarser mesh within a small, finite topological neighborhood, leading 

-18- 
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to very efficient implementations. Using a suitable choice of weights, such subdivision 
schemes can be designed to produce a smooth surface in the limit. Subdivision methods 
which result in limit surfaces whose curvature tensor is square integrable are especially 
appealing for geometrical modeling applications and for the purpose of thin-shell analysis. 

Many subdivision schemes have been proposed and studied extensively in the 
mathematical geometric modeling literature. The methods can be separated into two groups. 

Interpolating Schemes: The nodal positions of the coarser mesh are fixed, while 
only the nodal positions of new vertices are computed when going from a coarser to a finer 
mesh. Consequently, the nodal positions of the initial mesh, as well as any nodes produced 
during subdivision, interpolate the limit surface. Interpolating schemes exist for quadrilateral 
meshes and for triangular meshes. In both cases, the limit surfaces are but their curvatiu-es 
do not exist. Therefore, these schemes are not suitable for the special application of thin- shell 
analysis. 

Approximating Schemes: These schemes compute both new nodal positions for the 
newly created vertices, as well as for the vertices inherited from the coarser mesh, i.e., those 
which already carried nodal positions. Consequently, the nodal positions of the initial mesh 
are not samples of the final surface. The schemes of Catmull-Clark (E. CatmuU and J. Clark, 
"Recursively generated b- spline surfaces on arbitrary topological meshes," Comput. Aided 
Design, 10(6):350-355, 1978) and Doo-Sabin (D. Doo and M. Sabin, "Behaviour of recursive 
division surfaces near extraordinary points," Comput. Aided Design, 10(6): 356-360, 1978) 
fall into this class and operate on quadrilateral meshes. An approximating scheme for 
triangular meshes introduced by Loop (C. Loop, "Smooth subdivision surfaces based on 
triangles," Master's thesis. University of Utah, Department of Mathematics, 1987) is used in 
the present work. Loop's scheme produces limit surfaces which are globally except at a 
number of isolated points where they are only C\ However, their principal curvatures are 
square integrable, making them good candidates for thin-shell analysis. 

FIG. 3a shows a simple one-dimensional example of an approximating scheme 30. 
FIG. 3b shows an interpolating scheme 32. Both figures graphically show two subdivision 
steps based on an initial polygon. We assume that the polygon with the nodes x) is the result 
of subdivision step L In the subsequent refinement ^+1 for the approximating scheme 30, a 
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new vertex gets a nodal position jc*/^, which is the average of its two neighboring nodes jc* 



Since throughout the subdivision process all nodal positions are recomputed as weighted 
averages of nearby vertices, the resulting scheme does not interpolate the nodal positions of 
the control mesh, but rather approximates them. 

In contrast, the interpolating scheme 32 does not modify the coordinates of nodes 
existing at the previous refinement level: 



Repeating this process ad infinitum leads to smooth curves. In the case of the 
approximating scheme 30, these curves are actually cubic splines which are continuous. 
On the other hand, the interpolating scheme 32 is known as the 4pt scheme, and it can be 
shown that the resulting curves are C^'^ continuous. Since the entire subdivision process is 
linear, the resulting limit curves (or surfaces) are linear combinations of basis functions 
(sometimes referred to as "fundamental solutions" of the subdivision process). The compact 
support of the subdivision rules ensures that the basis functions are compactly supported as 
well. For the approximating scheme 30, this support extends two vertices to the left and two 
vertices to the right (and hence is a 2-ring), while the 4pt scheme has a support covering three 
vertices to each of the left and right (and hence is a 3-ring). Both schemes belong to the class 
of regular subdivision schemes since the weights are the same for every vertex and every 





The nodal positions of the existing vertices are recomputed as 




d^2I — •*'/ 



Only the newly generated vertices receive new nodal coordinates: 
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level. In the two-dimensional case, we will see that the weights will depend on the valence, 
i.e., the number of edges attached to a vertex, but are otherwise the same from level to level. 
These rules are also referred to as semi-regular. Similarly one can adapt the rules to non- 
smooth features such as boimdaries or creases, as well as other boundary conditions, using 
5 known techniques. 

In the one-dimensional case, subdivision schemes do not offer an important advantage 
since curves of the desired smoothness are easy to construct with traditional approaches such 
as Hermitian interpolation. In the two-dimensional, arbitrary-topology manifold setting, by 
contrast, subdivision methods offer significant advantages over other methods of smooth- 
10 surface construction. In fact, their original invention was motivated by the difficulties of 

constructing smooth-surface models of arbitrary-topology. For example, it is well known that 
arbitrary-topology surfaces built with traditional patches require up to sixth-order 

Q 

polynomials, leading to cumbersome computations and difficult-to-manage cross-patch 
^ continuity conditions. Additionally, the resulting degrees of freedom often lack physical 

M 15 meaning. However, in the subdivision setting, the only degrees of freedom are the nodal 
=J positions and the resulting surfaces are guaranteed to be smooth without the need to enforce 

^ cross-patch continuity conditions. 

'■5 

O In the following sections we introduce the refinement rules used in Loop's 

1=1] subdivision scheme for surfaces, and, based on the concept of the subdivision matrix, briefly 

20 discuss the basic ideas behind the scheme's smoothness analysis. This discussion will serve 

g ' ^ 

O to prepare the way for the parameterization of subdivision surfaces needed for the evaluation 

of positions, tangents, and curvatures. In particular, it is possible to evaluate these quantities 
exactly through the use of eigenanalysis without going to the limit. While we focus on 
Loop's scheme, the basic ideas and machinery apply equally well to other subdivision 
25 schemes. In particular, the CatmuU-Clark scheme, with quadrilateral elements, is a promising 
alternative for finite-element computations, since quadrilateral elements generally perform 
relatively better for very coarse meshes. 
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Refinement Rules 

In Loop's subdivision scheme, the control mesh and all refined meshes consist of 
triangles only. These are preferably refined by quadrisection (see FIG. 2). After the 
refinement step, the nodal positions of the refined mesh are computed as weighted averages 
of the nodal positions of the unrefined mesh. We distinguish two cases: new vertices 
associated with the edges of the coarser mesh, and old vertices of the coarse mesh. 

The coordinates of the newly generated nodes x\, x\, x\, on the edges of the 
previous mesh are computed as: 

^ 8 

whereby index / is to be understood in modulo arithmetic. The old vertices get new nodal 
positions according to: 

ajj+^ = (1 - Nw)xq + wx^ H + wx% 

where jc* are the nodal positions of the mesh at level k and jc*^* are the respective positions for 
the mesh ^+1 . The valence of the vertex, i.e., the number of edges incident on it, is denoted 
by N. Note that all newly generated vertices have valence 6 while only the vertices of the 
original mesh may have valence other than 6. We will refer to the former case (valence = 6) 
as regular and to the latter case as irregular. Equations 51 and 52 are depicted in symbolic 
form in FIGS. 4a and 4b as so-called subdivision masks 40, 41 . FIG. 4a shows a subdivision 

mask 40 for computing as a function of a parameter w. FIG. 4b shows a subdivision mask 

41 for computing Xj , ... jc] . 

At this stage it is not obvious how to choose the parameter w for the subdivision mask 
40 to get continuous surfaces. In the original scheme Loop proposed 



5 _ 1 
8 U 4 



27r\' 



As it turns out, other values for w also give smooth surfaces. For example. 
Warren's (J. Warren, "Subdivision methods for geometric design," Unpublished manuscript. 
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Department of Computer Science, Rice University, November 1995) choice for w is simpler 
to evaluate than that of Equation 53: 

3 3 
t£; = — for N>3 ^ = ^ ^ ^ (54) 

Although the choice for the weights appears somewhat arbitrary, the motivation for 
this choice will be discussed in the next section. In any case, the weights used by the 
subdivision scheme depend only on the connectivity of the mesh and are independent of the 
nodal positions. 

So far we have ignored the boundary of the mesh, where the subdivision rules need to 
be modified. Two choices are possible here. The first method considers the boundary as a 
one-dimensional curve and applies the rules of Eqs. 47 and 48 to any vertices on the 
boundary. Another method for the treatment of boundaries, which we employ in our 
implementation, was proposed by Schweitzer (J.E Schweitzer, Analysis and Application of 
Subdivision Surfaces, Ph.D. dissertation. Department of Computer Science and Engineering, 
University of Washington, Seattle, 1996). For each boundary edge, one temporary or 
artificial vertex is defined, after which the ordinary rules are applied. FIG. 5 gives a graphical 
example of a boundary approximation using temporary vertices 50, 52. The nodal positions 
of the temporary vertices 50, 52 are set to: 

This particular choice of temporary vertices effectively reproduces the one-dimensional 
subdivision rules. 

This simplified boundary treatment with temporary vertices leads to approximation 
errors for comers and boundary vertices with valence other than 4. However, these errors do 
not inhibit the convergence of the finite-element method. For modeling surfaces with creases 
or shell intersections, the edge rules can also be applied within the mesh on appropriately 
tagged edges. FIG. 6 shows an example of a geometrically more complex shape (a distributor 
cap) using such crease rules, demonstrating the ability of subdivision surfaces to model 
intricate shapes in a straightforward fashion. The entire shape is a single, piecewise smooth 
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surface. In this example, the original control mesh 60 is chosen to be least-squares optimal in 
terms of geometric approximation fidelity. FIG. 6 also shows a subdivided mesh 62 and the 
limit surface 64.. 

Convergence and Smoothness 

Subdivision methods for arbitrary-topology surfaces were introduced more than 20 
years ago, but not widely used in applications until the early 1990s when theoretical 
breakthroughs put their convergence and smoothness analysis on a solid mathematical 
foundation. If all vertices are regular (valence 6 for triangles or valence 4 for rectangular 
schemes), powerful Fourier methods can be used to establish smoothness properties. 
However, these techniques do not apply in the arbitrary-topology manifold setting, in which 
irregular vertices cannot be avoided. The critical question in this setting concerns the 
smoothness properties around irregular vertices. The key tool in this analysis is the local 
subdivision matrix and its structure. The fact that the analysis of these schemes can be 
reduced to the analysis of a local matrix is due to the local support of the basis functions. In 
other words, the behavior of the surface in a neighborhood of a node depends only on those 
basis functions whose support overlaps a neighborhood of the node. 

For Loop's scheme, the relevant neighborhood around a vertex is a 2-ring, i.e., all 
vertices which can be reached by traversing no more than two edges. The subdivision matrix 
expresses the linear relationship between the nodal positions in a 2-ring around a given vertex 
at level k and the nodal positions around the same vertex in the 2-ring at level k-^l . This 2- 
ring analysis is required to establish smoothness. Details of this approach, including 
necessary and sufficient conditions, can be found in the literature. Once the analytic 
properties have been established, quantities at a point, such as position or tangents, can be 
computed using an even smaller subdivision matrix which relates 1 -rings at level k to those at 
level k+l. In the following we will discuss only this simpler setting since it is the one which 
is needed for actual computations. 

Let A* be the vector of nodal positions of a vertex with valence A/^ and its 1-ring 

neighbors at level k, = ( jCq , jcf . . , jc^ ). Note that in the surface case we treat this vector 
as an A^+1 vector of three-dimensional vectors, while in the functional setting it would be an 

-24- 



Docket No.: 0661g/5^|pl/CIT3061 ' 

A^+l vector of scalars. We can now express the linear relationship between the nodal 
positions at level k and ^+1 with an (//+l)x(^+i) matrix 5: 

The entries of the matrix S are given by the subdivision rules (Eqs. 52 and 51). The 
study of the limit surface then amounts to examining 

- lim S^X^ 

From this it immediately follows that the limit surface cannot exist if S has an 
eigenvalue of modulus larger than 1 . Furthermore, it can be shown that it must have a single 
eigenvalue 1 and all other eigenvalues must have modulus strictly smaller than 1 . This 
property also implies that subdivision schemes are affinely invariant, /.e., an affine 
transformation applied to the nodal positions of the control mesh results in a limit surface 
having undergone the same affine transformation, which is a desirable property in practical 
applications. The associated right eigenvector is easily seen to be the vector of all I's. 

In the following we summarize the main results as needed for our finite-element 
setting. Assume that the subdivision matrix has a complete set of eigenvectors (this property 
holds for all schemes of practical interest, although it is not necessary for the analysis). Since 
the subdivision matrix S is not self-adjoint, let Rj and Lj be the right and left eigenvectors of 
S respectively, 

Lj • Ri = S jj 

and let be the associated eigenvalues in non-increasing magnitude order, with A-o = 1 . Using 
the eigenvalue decomposition, the nodal vector can be written as 

/=0 

with a°; = LfX^. Choosing this basis, Eq. 57 takes the simple form: 
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N iV 

- lim y a^jRj = Urn T X^a^jRf (60) 

From this representation and the facts that ^ = 1 and < 1 for /= 1,...^, it follows 
immediately that the center nodal position as well as all nodal positions in the 1-ring 
converge to 

X^^alRo (61) 

Since has the form ifg = (1, 1, 1) all the nodal positions in the 1-ring neighborhood 
actually converge to a single position: 

al = L^^ (62) 
For the original. Loop subdivision scheme, Lq is given as: 

r. = (i-;vM,.-,0 with i = (63) 

In practical terms this means that we can evaluate the limit position of the surface 
given at any finite subdivision level by simply applying the position mask 70 shown in FIG. 
7a. We can continue this analysis and compute tangent vectors in a similar way. Assume that 
the control nodal positions have been translated by so that the limit position for our 
selected vertex is the origin. The leading term of Eq. 60 is now controlled by the 
subdominant eigenvalues. Here we have X, = X,2 > ^3. which results in: 

lim -r- = a?Ri + alR2 (64) 

indicating that the nodal positions in the 1-ring all converge to a common plane spanned by 
the vectors a^^ and a**2. While these vectors are not necessarily orthogonal, they do span a 
plane for almost all initial configurations. An explicit formula for two tangent vectors to the 
limit surface is: 

ti = Li • X^ t2 = L2- X^ (65) 
whence the shell director follows as: 
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ti X t2 

^3 



1*1 X t2| 

In certain settings we may have X| > > Th^ above equations for tangent vectors 
continue to hold, requiring an argument only slightly more involved than the one above. For 
Loop's scheme, given by 

Li = (0. Ci, C2r • • 5 Cn) L2 = (0, 5i, 52, • • • , Sm) 

with 

27r(/-l) . 27r(/-l) 
C/ = cos St = sm — ^ . 

Again we find that a very simple computation allows us to compute the limit surface 
normal given the subdivision mesh at any level k. FIGS. 7b and 7c show two tangent masks 
72, 74 corresponding to Eqs. 65 and 67. 

So far we have only discussed the evaluation of position and tangents of the limit 
surface at parametric locations corresponding to control vertices. For numerical evaluation of 
the Kirchhoff-Love energy functional we need to evaluate these quantities - and curvature 
quantities - at quadrature points. This is particularly simple in the case of Loop's scheme 
(and similarly so for Catmull-Clark's scheme) since Loop's scheme actually generalizes 
quartic box splines (Catmull-Clark's scheme generalizes bi-cubic splines). If the valencies of 
the vertices of a given triangle are all equal to 6, the resulting piece of limit surface is exactly 
described by a single quartic box-spiine patch, for which very efficient evaluation schemes 
exist at arbitrary parameter locations. We call such a patch regular, FIG. 8 is a graphical 
example of a regular box-spline patch with 12 control points. Regular patches are controlled 
by 12 basis, or shape, functions (see Appendix A.0.1), since only their support overlaps the 
given patch. If a triangle is irregular, i.e., one of its vertices has valence other than 6 the 
resulting patch is not a quartic box spline. Arbitrary parameter locations may nevertheless be 
treated simply by the method described in the next section. 



Function Evaluation for Arbitrary Parameter Values 
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In this section we discuss the evaluation of function values and derivatives for Loop 
subdivision surfaces at arbitrary parameter locations. These function evaluations arise during 
the computation of element stiffness and force arrays by numerical quadrature, Eq. 46. 
Despite early attempts, the proper parameterization of subdivision surfaces has been until 
recently an unsolved problem in the vicinity of irregular vertices. Stam (J. Stam, "Fast 
evaluation of catmuU-clark subdivision surfaces at arbitrary parameter values," in Computer 
Graphics (SIGGRAPH V8 Proceedings), 1998; J. Stam, "Fast evaluation of loop triangular 
subdivision surfaces at arbitrary parameter values," in Computer Graphics (SIGGRAPH '98 
Proceedings, CD-ROM supplement) , 1998) presented an elegant solution to this problem 
based on the eigendecomposition of the subdivision matrix. 

A convenient local parameterization of the limit surface may be obtained as follows. 
For each triangle in the control mesh we choose {Q\ 9^) as two of its barycentric coordinates 
within their natural range 

T={{e\e% s. t. 5" € [0, 1]} (68) 

The triangle T in the (0\ 0^)-plane may be regarded as a master or standard element domain. 
It should be emphasized that this parameterization is defined locally for each element in the 
mesh. The entire discussion of parameterization and function evaluation may therefore be 
couched in local terms. 

Consequently we proceed to consider a generic element in the mesh and introduce a 
local numbering of the nodes lying in its immediate 1 -neighborhood. For regular patches 
such as depicted in FIG. 8, Loop's scheme leads to classical quartic box splines. Therefore, 
the local parameterization of the limit surface may be expressed in terms of box-spline shape 
functions, with the result: 

12 



where now the label / refers to the local numbering of the nodes. The precise form of the 
shape functions (9*, 0^) is given in Appendix A.0.1 . The surface within the shaded triangle 
in FIG. 8 is controlled by the 12 local control vertices. In contrast to Hermitian interpolation. 
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the surface is solely controlled by the position of these control vertices, and first- or second- 
order derivatives at the nodes are not utilized. The image of the standard domain T under the 
embedding (Eq. 69) constitutes the geometric domain of the element under consideration 
within the limit surface. The embedding may thus be regarded as a conventional 
isoparametric mapping from the standard domain Tonto Q, with {Q\ Q') playing the role of 
natural coordinates. 

The parameterization in Equation 69 may also be used for the displacement field. 
Locally one then has 

Uf,{e\,e') = f:N'{e\e')uj 

For fimction evaluation on irregular patches the mesh has to be subdivided until the 
parameter value of interest is interior to a regular patch. At that point the regular box-spline 
parameterization applies once again. It should be noted that the refinement is performed for 
parameter evaluation only. For simplicity we assume that irregular patches have one irregular 
vertex only. This restriction can always be met for arbitrary initial meshes through one step 
of subdivision, which has the effect of separating all irregular vertices. 

FIG. 9 illustrates the basic idea for function evaluation at arbitrary parameter values. 
It shows a patch 90 having a single vertex 91 of valence 7. After one level of subdivision this 
irregular patch 90 is divided into three regular sub-patches 93-95, and one irregular sub-patch 
96. If the desired parameter value lies within the regular sub-patches we may immediately 
evaluate the surface by using the canonical regular-patch evaluation routine. If our desired 
parameter location lies within the irregular sub-patch 96, we must repeatedly subdivide until 
desired parameter location falls within a regular patch. This can be achieved for any 
parameter value away from the irregular vertex; for the irregular vertex itself we have already 
described the evaluation procedure above. 

The action of the subdivision operator for the entire neighborhood defined by the sub- 
patches 93-96 can again be described by a matrix: 
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- AX"" 

The matrix y4 now has dimension (A/-I-12, and its entries can be derived from the 
subdivision rules as presented above. For the proposed shell element with one-point 
quadrature at the barycenter of the element, a single subdivision step is sufficient, since the 
sampling point lies in sub-patch 94. We define 12 selection vectors i*;, 1,...,12 of 
dimension (A^-f 12) which extract the 12 box-spline control points for sub-patch 94 from the 
iV+12 points of the refined mesh. The entries of Py are zeros and one depending on the indices 
of the initial and refined meshes. To evaluate the fiinction values in the three triangles 93-95 
with the box-spline shape ftmctions N^, a coordinate transformation must be performed. The 

relation between the coordinates (0 \ 9 ^ ) of the original triangles and the coordinates 

(9 \ ) of the refined triangles can be established from the refinement pattern shown in FIG. 
10. For the center of sub-patch 94 we have the following relation: 

The function values and derivatives for sub-patch 94 can now be evaluated using the 
interpolation rule: 

12 

x{9\J^) = J^N^{9\9^)PjAX'' 

7=1 

Differentiation of the above equation yields the partial derivatives of x at the desired 
location. In these calculations, the jacobian matrix of the simple mapping (Eq. 73) between 

the coordinates (0\0^) and (0^0^) has to be included. The result is: 

1=1 

The jacobian of the embedding required in the quadrature rule (Eq. 46) can in turn be 
computed directly from through Eq. 8. Second-order derivatives such as required for the 
computation of bending strains follow simply by direct differentiation of the interpolation 
rule (Eq. 74). 
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The evaluation procedure just described is a simplified version of the scheme 
presented by Stam, cited above, for function evaluation at arbitrary parameter values. Since 
we are only interested in a one-point quadrature rule, a single subdivision step is sufficient. 
For other rules one may need to subdivide further. Stam describes the general case in which 
5 the eigendecomposition of^ is exploited to perform the implied subdivision steps through 
suitable powers of A in the eigenbasis, a simple and efficient procedure. 

Boundary Conditions 

Due to the local support of the subdivision scheme, the boundary displacements are 
influenced only by the nodal positions in the 1-neighboorhood of the boundary, /.e., the first 
10 layer of vertices inside the domain as well as a collection of artificial or "ghosf vertices just 
outside the domain. For example, for built-in boundary conditions the displacements of the 
nodes inside, outside, and on the boundary must be zero. The deformation of the boundary 
computed v^th the limit formula as described above shows that this results in zero 
f=j displacements and in fixed tangents. Other boundary conditions can be accommodated in a 

y 15 similar way. FIGS. 1 la, 1 lb, and 1 Ic show various boundary conditions for the cases of 
y displacements and rotations fixed (FIG. 1 la), displacement and rotations free (FIG. 1 lb), and 

;L, displacements fixed and rotations free (FIG. 1 Ic). 

N= The enforcement of prescribed boundary displacements is equally straightforward. 

ru 

•ff\ General displacement boundary conditions may be formulated with the aid of a local 

^ 20 reference frame defined by the surface normal and the tangent to the boundary of the shell. 

The prescribed displacement boundary conditions may be regarded as linear 
constraints on the displacement field and treated accordingly during the solution process. 
Such linear constraints may be enforced by a variety of known methods. In all the 
calculations discussed in the examples presented here, the displacement boundary conditions 
25 are enforced by the penalty method, with a penalty stiffness equal to 1 00 times the maximum 
diagonal component of the stiffness matrix. 

Computer Implementation of Subdivision-Based Shell Elements 
The current implementation of subdivision-based shell elements in accordance with 
the invention requires a few data structures in addition to those typical of conventional finite 
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element methods. For example, tables of element adjacencies are used for the computation of 
element arrays. In our particular implementation, we use a Triangle C-structure which 
stores the following information for each element of the mesh: 

typedef struct __T { 
5 /* pointer to the element vertices */ 

/* pointer to the neighbor triangles */ 

/* array of pointers to the vertices in the 1- neighborhood */ 
/* number of nodes in the 1 -neighborhood */ 

/* information about the boundary conditions of the edges */ 
10 } Triangle; 

The preferred method for performing finite element analysis using subdivision basis 
functions can be sunmiarized functionally as follows, referring to FIG. 12: 

1 . Input a mesh comprising a set of data points each having connectivity to neighboring 
data points, the mesh defining the physical state of an object or system (STEP 120); 

O 

15 2. Specify an initial state for the mesh (STEP 122); 

."k 3. Define a set of linear differential equations comprising a stiffness matrix (which 

represents the material and structural properties of the object or system) and an 

ft it 

Ln external forcing vector (which defines load conditions), at least one such equation 

" having a fourth order differential operator (this last requirement ensures that second 

H 20 order derivatives must exist, a requirement for functions) (STEP 124); 

fy 4. Solve the set of linear equations as applied to the mesh (STEP 126); and 

01 

Q 5. Output the solution to the set of linear equations as defining a modification of the 

~ initial state of the mesh based on the stiffness matrix and in response to the external 

forcing vector (STEP 128). 

25 This procedure can be implemented by programming the relevant equations set forth 

above in light of known techniques for finite element analysis. In more detail, referring to 
FIG. 13, one analytical solution process consists of the foUov^ng steps: 

1 . Model the geometry of a thin shell vAth a subdivision control mesh (STEP 1 30); 

2. Subdivide the initial subdivision control mesh once in order to separate the irregular 
30 vertices (the initial mesh could include triangles v^th more than one irregular vertex - 

in order to limit the number of possible mesh patterns during the parameter 
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evaluation, we refine the mesh once so that each triangle has at most one irregular 
vertex) (STEP 131); 

3. Generate artificial nodes and elements at the boundaries of the mesh (STEP 132); 

4. Find the 1 -neighborhood of each element of the mesh and gather the local nodes in 
5 accordance with the local numbering convention. (STEP 133); 

5. Create local coordinate systems at the nodes if necessary (for the creation of the local 
coordinate systems, this embodiment employs the limit formula of Eq. 66 for the shell 
normal) (STEP 134); 

6. Compute and assemble the element stiffiiess matrices and force arrays (STEP 135); 
10 7. Introduce displacement boundary condition constraints (STEP 136); 

8. Solve the system of equations (STEP 137); 

_ 9. Compute the limit positions of the nodes (Eq. 61) (STEP 138); and 

B 

10. Output a graphic or textual description of the deformed geometry of the modeled shell 
based on the computed limit positions of the nodes (see examples below for graphical 
15 output of deformed shell geometries) (STEP 139). 

y The elements under consideration here have three nodes and three displacement 

Q degrees of freedom per node. In a preferably embodiment, we use one-point quadrature with 

II, the sole quadrature point in the rule located at the barycenter of the element, i.e., at G^^ ==1/3 

'iU 

and 0^ =1/3. The corresponding weight \sw= 1/2, which, evidently, is the area of the 
20 standard triangle T. 

In sum, the preceding developments lead to the definition of a class of fiilly 
conforming "C*" triangular elements containing three nodes and one quadrature point. This 
combination of attributes, namely, the low order of interpolation and quadrature required, 
render the element particularly attractive from the standpoint of computational efficiency. As 
25 mentioned earlier, subdivision surfaces may also be used to define four-node square shell 
elements. The finite element analysis preferably uses subdivision basis functions as shape 
functions, but may in general use any suitably smooth shape ftmction. 

The invention may be implemented in hardware or software, or a combination of both 
(e.g., programmable logic arrays). Unless otherwise specified, the algorithms included as part 
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of the invention are not inherently related to any particular computer or other apparatus. In 
particular, various general purpose machines may be used with programs written in 
accordance with the teachings herein, or it may be more convenient to construct more 
specialized apparatus to perform the required method steps. However, preferably, the 
invention is implemented in one or more computer programs executing on programmable 
systems each comprising at least one processor, at least one data storage system (including 
volatile and non-volatile memory and/or storage elements), at least one input device, and at 
least one output device. The program code is executed on the processors to perform the 
functions described herein. 

Each such program may be implemented in any desired computer language (including 
machine, assembly, or high level procedural, logical, or object oriented programming 
languages) to communicate with a computer system. In any case, the language may be a 
compiled or interpreted language. 

Each such computer program is preferably stored on a storage media or device (e.g., 
solid state, magnetic, or optical media) readable by a general or special purpose 
programmable computer, for configuring and operating the computer when the storage media 
or device is read by the computer to perform the procedures described herein. The inventive 
system may also be considered to be implemented as a computer-readable storage medium, 
configured with a computer program, where the storage medium so configured causes a 
computer to operate in a specific and predefined manner to perform the functions described 
herein. 

Examples 

We proceed to establish the convergence characteristics of the inventive method by 
running the obstacle course of test cases proposed by Belytschko et al (7! Belytschko, H. 
Stolarski, W.K. Liu, N. Carpenter, and J.S.J. Ong, "Stress projection for membrane and shear 
locking in shell finite-elements," Comput. Methods Appl Meek Engrg,, 51:221-258, 1985). 
The shells in these tests cases take simple spherical or cylindrical shapes which can be 
readily described analytically. A preliminary step in the calculations is to approximate the 
exact surface of the shell by a subdivision surface in accordance with the invention. Several 
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methods of approximation are possible, including (1) least-squares approximation of the 
exact surface by the limit surface, and (2) placement of the control-mesh nodes on the exact 
surface. 

However, theoretical considerations and our own numerical tests show that the error 
5 incurred in the approximation of the shell geometry is of higher order than the finite-element 
error, and consequently both methods of approximation result in the same convergence rates 
asymptotically. It should also be noted that the question of geometrical approximation is 
rendered moot within an integrated computer-aided geometrical design (CAGD)-finite- 
element analysis framework. In this environment, the subdivision surface generated by the 
10 CAGD module becomes the true shell surface to be analyzed by the finite-element analysis 
module. 

We note for fiirther reference that a strictly C' finite-element method for Kirchhoff- 
Love shell theory containing the complete set of third-order polynomials within its 
interpolation satisfies the error bound: 

C 

luH-u\2.Si < -j^r\^kn (75) 

15 where u is the exact solution, u^^ is the finite-element solution, C is a constant, N is the 

number of nodes in the mesh, || ^ is the standard norm over H\Q; R% or "energy norm," 
1 13 Q is the standard semi-norm over //^(Q; 7?^). 

The central question to be ascertained now is whether the method developed in the 
foregoing exhibits the optimal convergence rate implied by the bound (Eq. 75). All the 
20 calculations described subsequently are carried out with one-point quadrature. The successive 
mesh refinements considered in convergence studies are obtained by regular refinement. 

We additionally compare the performance of the proposed approach with that of two 
other shell elements: 

ASM: The assumed-strain four-node shell element of Simo et ai (J.C. Simo, D.D. 
25 Fox, and M.S. Rifai, "On a stress resultant geometrically exact shell model, part ii: The linear 
theory; computational aspects," Comput. Methods Appl Meek Engrg., 73:53-92, 1989). 
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DKT-CST: A flat 3-node element with no membrane/bending coupling. In this 
element the discrete Kirchhoff triangle (DKT) formulation of Batoz (J.L. Batoz, "An explicit 
formulation for an efficient triangular plate-bending element". Internal J, Numer. Methods 
Engrg., 18:1077-1089, 1982) is utilized for the bending response, and standard constant 
strain interpolation is used for the membrane response. The resulting element has six degrees 
of freedom per node. 

Rectangular Plate 

As a first example, we consider the simple case of a square plate under uniform 
loading /? = 1 . FIG. 14 shows an example of such a rectangular plate 140, where the length of 
the plate 140 is Z = 100 and the thickness is = 1 . These dimensions place the plate well 
within the scope of KirchhofTs theory. The Young's modulus \sE= W and the Poisson's 
ratio is V = 0. In order to test the treatment of displacement boundary conditions described 
above, we analyze both a simply-supported and a clamped plate. The entire plate is 
discretized into finite elements with no account taken of the synmietry of the plate. A typical 
mesh 142 used in the calculations is shown in FIG. 14. The artificial or ghost nodes used to 
enforce the boundary conditions are not shown in the figure. 

FIGS. 15a and 15b respectively show the computed limit surfaces of the simply- 
supported plate 150 and the clamped plate 152 following deformation (deflections are scaled 
differently in both cases). The high degree of smoothness of these deflected shapes is 
noteworthy. An appealing feature of the problem under consideration is that it is amenable to 
an exact analytical solution. For instance, the maximum displacement at the center of the 
plate is found to be: u^^ ^ 0.487 for the simply-supported plate 150, and u^^ « 0.151 for the 
clamped plate 152. It therefore follows that the solution error can be computed exactly for 
this problem. FIGS. 16a and 16b respectively show graphs of the convergence of the 
computed maximum-displacement 160 and energy-norm 162 errors as a function of the 
number of degrees of freedom. In all cases, the optimal convergence rate 0{N^) is attained in 
the energy norm, which attests to the good convergence properties of the method. These 
results, and those presented subsequently, also demonstrate the sufficiency of the one-point 
quadrature rule. 
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Scordelis-Lo Roof 

The Scordelis-Lo Roof is a membrane-stress dominated problem and, as such, it 
provides a useful test of the ability of the finite-element element method to represent complex 
states of membrane strain. The problem concems an open cylindrical shell loaded by gravity 
forces. FIG. 17 shows an example of a Scordelis-Lo Roof cylindrical shell 170. In our 
calculations, the length of the cylinder is I = 50; its radius is 7? = 25; the angle subtended by 
the roof \s^ = 80°; the thickness is /z = 0.25; the Young's modulus is £ = 4.32 iO^ ; and the 
Poisson's ratio is v = 0. 

FIGS. 18a and 18b respectively show the undeformed control mesh 180 and the 
deformed limit surface 182. FIG. 19 shows a convergence plot for the maximum 
displacement of the deformed limit surface 182. The displacements are normalized by the 
value 0.3024 as given by Belytschko, cited above. The excellent convergence characteristics 
of the method are evident fi*om the figure. In particular, the subdivision element outperforms 
both the assumed-strain and the DKT elements. 

Pinched Cylinder 

The pinched-cylinder problem tests the inventive method's ability to deal with 
inextensional bending modes and complex membrane states. The problem concems a 
cylindrical shell pinched under the action of two diametrically opposite unit loads located 
within the middle section of the shell. We consider two cases: free-end boundary conditions; 
and ends constrained by two rigid diaphragms. FIG. 20a shows an example of a cylindrical 
shell 200. The length of the cylinder 200 is L = 600; the radius is R = 300; the thickness is h 
= 3; the Young's modulus is E = 3x 10^; and the Poisson's ratio is v = 0.3. 

A control-mesh node is placed at the point of application of the loads. It is interesting 
to note, however, that, owing to the nonlocal character of the shape functions, the point load 
is spread over several nodes. This is in contrast to other methods, e.g., those based on 
Hermitian interpolation. It should be carefully noted that the shape functions developed here 
possess the requisite partition-of-unity property and, in consequence, the resultant of all nodal 
forces exactly matches the applied load. 
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FIG. 20b shows a typical control mesh 202 for the cylinder of FIG. 20a. FIG. 20c 
shows the computed deformed limit surface 204, with the load point 206 on one side. Here 
again, the high degree of smoothness of the deformed surface, which attests to the 
conforming nature of the method, is noteworthy. FIGS. 21a and 21b respectively show 
displacement-convergence results for the rigid diaphragm case and the free-end case. In the 
rigid diaphragm case, we monitor the displacements under the loads, which are normalized 
by the analytical solution of 1 .82488 x 10 The excellent convergence properties of the 
method are evident from the figure. As may be seen, the subdivision element converges faster 
than both the assumed-strain and the DKT elements. The analytical solution for the free-end 
case has been given by Timoshenko (S. Timoshenko, Theory of Plates and Shells, 
McGraw-Hill Book Company Inc., New York-London, 1940). Timoshenko's solution gives a 
value of 4.520 xlO"^ for the displacements under the loads, and a value of 4.156 x 10^ for the 
change in diameter at the free ends. The convergence of this latter quantity is shown FIG. 
21b. Here again, the robust convergence characteristics of the numerical solution is 
noteworthy. 

Hemispherical Shell 

The case of a spherical shell provides an example of a surface which cannot be 
triangulated without the inclusion of irregular nodes, or nodes of valence different from 6. 
Under these conditions, a naive implementation of a box-spline-based finite-element method 
necessarily breaks down, which illustrates the need for the more general treatment developed 
here. We consider a shell of hemispherical shape deforming imder the action of four point 
loads acting on its edge. FIG. 22a shows an example of a hemispherical shell 220. The radius 
of the hemisphere is i? = 10; the thickness \sh = 0.04 ; the Young's modulus isE = 6.825 10^; 
and the Poisson's ratio is v = 0.3. The edge of the shell is free. The applied loads have a 
magnitude F = 2 and define two pairs of diametrically opposite loads alternating in sign at 
90°. 

FIGS. 22b and 22c respectively show a typical control mesh 222 and the 
corresponding deformed limit surface 224. The presence of irregular nodes in the control 
mesh 222 should be carefiilly noted. FIG. 23 shows the convergence of the radial normalized 
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displacement under the applied loads. The displacements are normalized by the exact 
solution 0.0924 as given by Belytschko, cited above. The hemispherical shell is a standard 
test for assessing an element's ability to represent inextensional deformations. Critical to the 
performance of the element is its ability to bend without developing parasitic membrane 
strains. Accordingly, the flat DKT element with uncoupled membrane and bending stiffiiess 
performs particularly well in the pinched-hemisphere shell problem. However, a careful 
examination of FIG. 23 reveals that the subdivision element attains convergence faster than 
even the DKT element. The excellent convergence characteristics of the method are again 
evident from the figure. 

Conclusions 

We have presented a new paradigm for thin-shell finite-element analysis based on the 
use of subdivision surfaces for (1) describing the geometry of the shell in its undeformed 
configuration, and (2) generating smooth interpolated displacement fields possessing 
bounded energy within the strict framework of the Kirchhoff-Love theory of thin shells. 
Several salient attributes of the inventive interpolation scheme bear emphasis: 

• The undeformed and deformed surfaces of the shell, or equivalently the displacement 
field thereof, follow by the recursive application of local subdivision rules to nodal 
data defined on a control triangulation of the surface. The particular subdivision 
strategy adopted here is Loop's scheme. Whether directly at regular nodes of valency 
6, or following an appropriate number of subdivision steps in the case of irregular 
nodes, the limit surface can be described locally by quartic box splines. Subdivision 
rules may also be defined for square meshes. The subdivision rules may also be 
generalized to account for creases and displacement boundary conditions. 

• The limit surfaces obtained by subdivision, and the displacement fields they support, 
are except at isolated extraordinary points. The displacement field is guaranteed to 
be ff^ i.e., to be square-summable and to have first and second square-summable 
derivatives. In consequence, the interpolated displacement fields have a finite 
Kirchhoff-Love energy. In the parlance of the finite-element method, the proposed 
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interpolation scheme is strictly C* and therefore meets the convergence requirements 
of the displacement finite-element method. 

• The triangles in the control mesh induce subregions on the limit surface which may 
be regarded as bona fide finite elements. In particular, the energy - and all other 
extensive properties - of the shell may be computed as a sum of integrals extended to 
the domains of the elements. These element integrals may conveniently be evaluated 
by numerical quadrature without compromising the order of convergence of the 
method. A one-point quadrature rule is sufficient for the computation of element 
integrals resulting fi-om Loop's subdivision scheme. 

• The displacement field of the shell is interpolated fi-om nodal displacements only. In 
particular, no nodal rotations are used in the interpolation. However, the interpolation 
scheme developed here differs from conventional C° finite-element interpolation in a 
crucial respect: the displacement field over an element depends on the displacements 
at the three nodes of the element and on the displacements of the 1-ring of nodes 
connected to the element. In this sense, the interpolation rule is nonlocaL However, 
the use of subdivision surfaces ensures that all the local displacement fields defined 
over overlapping patches combine conformingly to define one single limit surface. 

In sum, subdivision surfaces enable the finite-element analysis of thiri shells within 
the strict confines of Kirchhoff-Love theory while meeting all the convergence requirements 
of the displacement finite-element method. In particular the ability to couch the analysis 
within the framework of Kirchhoff-Love theory entirely sidesteps the difficulties associated 
with the use of methods in the limit of very thin shells. In a particularly pleasing way, 
subdivision surfaces enable the return to the most basic-and fimdamental-of finite element 
approaches, namely, constrained energy minimization over a suitable subspace of 
interpolated displacement fields, or Rayleigh-Ritz approximation. Finite-element methods 
formulated in accordance with this prescription satisfy the orthogonality property, i.e., the 
error function is orthogonal to the space of finite-element interpolants; and possess the best- 
approximation property, i.e., the energy norm of the error is minimized by the finite-element 
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solution. These properties render the basic finite-element method exceedingly robust and 
account for much of its success. Our numerical experiments show that the approach proposed 
here does indeed lead to the optimal convergence rate predicted by finite-element theory. 

Another key advantage afforded by the approach developed here is that subdivision 
surfaces provide a common representational paradigm for both solid modeling and shell 
analysis, with the attendant unification of traditionally heterogeneous software tools. By 
virtue of this unification, surface geometries generated by a computer-aided geometry design 
(CAGD) module can be directly utilized by the shell-analysis module without the need for 
any intervening geometrical manipulation. As a consequence, high-level algorithms 
developed in the field of computer aided geometric design can be integrated simply into the 
shell analysis software. 

A number of possible extensions of the theory and applications of the theory are 
worth mentioning. Firstly, recursive subdivision provides an effective basis for mesh 
adaptation. By retaining the hierarchy of finite-element representations generated by 
subdivision, the application of multiresolution methods and related techniques, such as 
wavelets, becomes straightforward. Indeed, the application of wavelet methods to the 
description of complex and intricate geometries has already been extensively pursued within 
the field of computer graphics. Finally, the extension of the proposed approach to the 
nonlinear range appears straightforward. In this particular context, the sole use of nodal 
displacements in the interpolation is expected to simplify the solution procedure by 
eliminating the need for introducing complex schemes for the nonsingular parametrization of 
the shell director. 

A number of embodiments of the present invention have been described. Neverthe- 
less, it will be understood that various modifications may be made without departing fi-om the 
spirit and scope of the invention. For example, while the discussion above has focused 
applied forces as the loading condition, the loading conditions may include thermal loading 
Accordingly, other embodiments are within the scope of the following claims. 
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Appendix 

A. 0. 1 Regular Patches 

For regular patches the shape functions are given by the 12 box-spline basis functions. 
The loczd numbering of the nodes adopted here is as in FIG. 8. 

N2 = ^iu^ + 2u^w) 

N4 = ^ (6w* + 24u^w + 2Au^w'^ + 8mw^ + 10" + 24u^y + QQu^vw + ZQuvv? 

+ Qvw^ + 24u^v^ + 36uv'^w + 12i;^w^ + Suv^ + 6v^w + u") 

N5 = ^ ^^^^ 12w^w;^ + 6mu^ + + 2u^v + Qu^vw + 6m?;w^ + 2?;u'^) 

Ne = ^i2uv' + v'') 

N7 = ^ (w* + 6it''w + 12tf^w;^ + 6Mxy^ + + 8w^t; + 36w^i;w + 36uvw^ + 8vw^ 

J. ^ 

+ 2Au\'^ + eOuv^w + 2Av^w^ + 2Auv^ + 247;^t/; + Gv'') 
ATg = — (u* + Sti'^w + 24'u^zi;^ + 24ww^ + Qw"^ + Gw^t; + ZQv^vw + GOwu'u;^ 

+ 24t;it;^ + l2uW + ZQuv^w + 24?;^m;^ + 6wi;^ + 8v^w + y^) 
iVg = ^[2uw^ + w'') 

N,o - ^(2^;2^' + i;^) 

A'u = ^{2uw^ + w'^ + 6uvw^ + 6vw^ + 6uv'^w + I2v'^w^ + 2uv''^ + 6v^w + v'^) 
N12 = ^{w' + 2vw^) (76) 
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where the bary centric coordinates (w, v, w) are subject to the constraint: 

u + V + w = 1 

The local curvilinear coordinates (9*, 6^) for the element may be identified with the 
barycentric coordinates (v, w). 

AM. 2 Irregular Patches 

As discussed above, a closed-form representation for the shape functions is not 
available at irregular vertices. The shell surface within each element is, however, completely 
described by the nodal positions of the element and its 1-ring. For the one-point quadrature 
rule used in calculations, the regular patch configuration is recovered after the application of 
one single subdivision step. In the following we give a more general version of the function 
evaluation scheme for arbitrary parameter values. The number of subdivision steps required 
and the attendant coordinate transformations can be computed by the following algorithm due 
to Stam (J. Stam, "Fast evaluation of loop triangular subdivision surfaces at arbitrary param- 
eter values," in Computer Graphics (SIGGRAPH '98 Proceedings, CD-ROM supplement), 
1998): 

/* determine the necessary number of subdivisions */ 
u = v+w; 

k = -loglO (u) /loglO (2 .0) ; 
k = ceil(k) ; 

pow2 = pow(2.0, (double) (k-l) ) / 

/* determine in which domain (v, w) lies */ 

V * = po w2 ; w * = pow2 ; 

if (v > 0.5) { 

V - 2.0*v-1.0; w = 2.0*w; 
whpa = 1; 

} else if (w > 0.5) ( 

V = 2.0*v; w = 2.0*w-1.0; 
whpa = 3; 

} else { 

V = 1.0 - 2.0*v; w - 1.0-2.0*w; 
whpa = 2; 

} 
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The function value at any parameter value is given by: 



NP 



7=1 

which generalizes Eq. 73. Here again, the vectors Pf extract the control nodes of the box- 
spline patch. The variable whpa in the above code gives the number of the subpatch, which 
contains the coordinates (v, w), and the variable k the power of the matrix y4. For Loop's 
scheme, the subdivision matrix^ has the following form: 



A = 



with 
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It should be noted that the dimensions of the matrices S^^ and depend on the 
valence of the irregular vertex. 
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A. 0,3 Membrane and Bending Strain Matrices 

The membrane and bending-strain matrices take the form; 



and 



A^^.i 0.1 • ei ai ■ 62 N^,i ai ■ 63 

N^,2 02 • Ci N',2 02 • 62 A^^2 a2 • €3 

(iV^2 ai + iV^l 02) • ei (A^^2 ai + yV^i 02) • 62 (A^',2 ai + AT^^ 02) 



B[-ei Bi'e2 B^-e^ 
Bi -ei Bi-e2 63 



respectively. In-the above expressions, (e„ tfjj^a) are the basis vectors of an orthonormal- 
coordinate reference frame, and 

+ 0,3 • ai^i{N^,i a2 X 03 + N^,2 03 x 03)] 



B2 = -A^^,22 az + -7=[A''^.i 02,2 X 02 + A^',2 ai x 02,2 

+ O3 • 02,2 (AT^,! O2 X 03 + N^,2 03 X 03)] 

= -^'^12 O3 + 4=[^'^^1 '*1.2 X 02 + Ar^2 Ol X Oi,2 

+ 03 • Oi,2(A/'^,i O2 X O3 + A/'^,2 03 X 03)] 
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